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Motivated by the example of a curved waveguide embedded in a photonic crystal, we examine 
the effects of geometry in a "quantum channel" of parabolic form. We study the linear case and 
derive exact as well as approximate expressions for the eigenvalues and eigenfunctions of the linear 
problem. We then proceed to the nonlinear setting and its stationary states in a number of limiting 
cases that allow for analytical treatment. The results of our analysis are used as initial conditions in 
direct numerical simulations of the nonlinear problem and localized excitations are found to persist, 
as well as to have interesting relaxational dynamics. Analogies of the present problem in contexts 
related to atomic physics and particularly to Bose- Einstein condensation are discussed. 

I. INTRODUCTION 

Localization phenomena are widely recognized as key to understanding the excitation dynamics in many physical 
contexts such as light propagation, charge and energy transport in condensed-matter physics and biophysics, and 
Bose-Einstein condensation of dilute atomic gases [1,2]. Recent advances in micro-structuring technology have made 
it possible to fabricate various low-dimensional systems with complicated geometry. Examples are photonic crystals 
with embedded defect structures such as microcavities, waveguides and waveguide bends [2-5], narrow structures 
(quantum dots and channels) formed at semiconductor heterostructures [6-8], magnetic nanodisks, dots and rings 
[9-11]. 

On the other hand, it is well known that the wave equation subject to Dirichlet boundary conditions has bound 
states in straight channels of variable width [12-14], and in curved channels of constant cross-section [15,16]. Spectral 
and transport characteristics of quantum electron channels [17] and waveguides in photonic crystal [3] are in essential 
ways modified by the existence of segments with finite curvature. The two-dimensional Laplacian operator supported 
by an infinite curve which is asymptotically straight has at least one bound state below the threshold of the continuum 
spectrum, as was recently proved in [18]. The appearance of an effective attractive potential in the wave equation 
is due to constraining quantum particles from higher to lower dimensional manifolds [19-21]. Curvature induced 
bound-state energies and and corresponding wave functions were studied in [22] . 

Until recently there have been few theoretical and numerical studies of the effect of curvature on properties of 
nonlinear excitations. Nonlinear whispering gallery modes for a nonlinear Maxwell equation in microdisks were 
investigated in [23] ; the excitation of whispering-gallery- type electromagnetic modes by a moving fluxon in an annular 
Josephson junction was shown in [24]. Nonlinear localized modes in two-dimensional photonic crystal waveguides 
were studied in [25]. A curved chain of nonlinear oscillators was considered in [26] and it was shown that the 
interplay of curvature and nonlinearity leads to a symmetry breaking when an asymmetric stationary state becomes 
energetically more favorable than a symmetric stationary state. Propagation of Bose-Einstein condensates in magnetic 
waveguides was experimentally demonstrated quite recently in [27]; single- mode propagation was observed along 
homogeneous segments of the waveguide, while geometric deformations of the microfabricated wires led to strong 
transverse excitations. 

Motivated by the experimental relevance of the above mentioned geometric deformations, in the present work we aim 
at investigating nonlinear excitations in a prototypical setup incorporating such phenomena. As our case example, we 
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will examine an infinitely narrow curved nonlinear waveguide (channel) embedded in two-dimensional linear medium 
(idle region). 

Our presentation will proceed as follows: in Section II, we set up the mathematical model of interest and examine its 
general properties and equations of motion. In Section III, we study the linear case and present its explicit solutions 
for the bound states, as well as for the corresponding eigenvalues. In section IV, we focus on the nonlinear problem, 
while in section V we supplement our analysis with numerical results. Finally, in section VI, we summarize our 
findings and present our conclusions. The appendix presents some additional technical details on the solution of the 
linear problem and its Green's function formulation. 



II. SYSTEM AND EQUATIONS OF MOTION 

Our model is described by the Hamiltonian 



H 



OO 

J {|V^| 2 - v (H 2 - A^\ 4 ) S(y - f(x))}dxdy, (1) 



where ip{r 7 t) is the complex amplitude function, f = (x,y), V 2 = d 2 + d 2 , v is the energy difference between the 
quantum channel and the passive region (refractive index difference in the case of photonic crystals and waveguides), 
the coefficient A characterizes the nonlinearity of the medium, e.g., the nonlinear corrections to the refractive index 
of the photonic band-gap materials, or self-interaction of the quasi-particles in the quantum channel. The function 
y = f(x) gives the shape of the channel. From the Hamiltonian, we obtain the equation of motion in the form 

id t i>(r, t) + V 2 ^ + v8{y- f(x))F(\^\ 2 ) V> = 0, (2) 

where the function -^(IV^I 2 ) is given by 

F=l-A\^\ 2 . (3) 

Equation (2) has as integrals of motion the Hamiltonian (1) and the (L 2 ) norm (referred to e.g., as the number of 
atoms in BEC or power in nonlinear optics) 



N 



OO 

f Wdxdy. (4) 



It is well known (see e.g., [28]) that if the edges of the Frenet trihedron (the tangent, the principal normal and the 
binormal) at a given point are considered as the axes of a Cartesian coordinate system, then the equation of the curve 
in a neighbourhood of this point has the form 

x = s + ---, y = -s 2 + ---, z = - — s J + • • • (5) 
2 o 

where s is the arclength, k is the curvature and r is the torsion of the curve at this point. Thus if the curvature of 
the plane curve is not too large one can represent it as a parabola 

K 2 



y=2 x - ( 6 ) 



In this case Eq. (2) takes the form 



z3^ + V 2 V> + ^%-!^)F(|Vf)V = 0, (7) 
where R — 1/k is the maximum radius of curvature of the curve. It is convenient to use the parabolic coordinates 

X =R> y =2 + 2R {u - v) - (8) 



2 



The coordinate lines are two orthogonal families of confocal parabolas, with axes along the y axis. These lines are 
given by 



2 2 

R% = 2y - R + v 2 /R, R^ = -2y + R + u 2 / R 



or 



u/VR = ±\y-^ + J(y-^) 2 + x 2 , 



v/^R=]J-y + § + \j(y-^Y + x 2 . (9) 

The variable u is allowed to range from — oo to oo, whereas v is positive. Introducing Eqs. (8) into Eq. (7) and using 
the properties of parabolic coordinates (see e.g., [29]), we obtain: 

* ^2 (" 2 + d ^ + ( d u + d v) 4> + v5(v-R) F(\^j\ 2 ) i> = Q. (10) 
Using the Fourier transform with respect to t, 

oo 

^=2^1 e^(«,M)*, (11) 

— OO 

where the bars denote the Fourier transformed quantities, one can represent Eq. (10) in the form 

-^{u 2 + v 2 ) j+(d 2 u + d 2 v ) 4, + v5(v-R) iW)V>K ») = o- (12) 
Equation (10) can, in turn, be expressed in the form of the integral equation 

oo oo 

${u,v,w) = v j du' j dv' G (u, v; u', v') 6(v' - R) F(\tp\ 2 ) ip(u', v', u), (13) 

-oo 

where the Green's function G (u, v; u', v') satisfies the equation 

(d 2 u + d 2 ) G -^(u 2 +v 2 )G = -S(u u') 5{v v'), (14) 
and has the form (see Appendix for details) 

, CO 

G(u,v;u',v') = Fn(u)F n (u') {V n (v) U n (v') 6(v' - v) + U n (v)V n (v') 6(v-v')}. (15) 

n=0 

Here 

F n {u) = a n e-^ u2 ' 2R H n (uuj 1 l i R- 1 l 2 ), n = 0,l,2,.., (16) 
where H n (z) is the Hermite polynomial [30], a n = [u 1 / 2 / Rtt) 1 ^ (2"n!) -1 ^ 2 is the normalization constant, and 

U n {v)=u{n+ l - 1 vV2uj 1 ' i R- 1 l 2 ^ 1 (17) 

with V(a, x) and U(a,x) being the Weber parabolic cylinder functions [30]. 

It can then be seen from Eqs. (13) and (15) that the Fourier transformed channel wave function 



<t>{u,w) = ifj(u,v,ui) 



v=R 



satisfies the equation 



u,u>) = vl>- 1 '*^^ J2 I F n (u)F n (u')V n (R)U n (R)F(\^)^u',u;)du\ 



n=0 



or equivalently the equation 

: , . oo 

n J 



n=0 



F n (u) F n (u') - , 



The wave function ip(u, v, ui) may be represented in terms of the channel wave function <fi(u, uj) as follows: 

rm 



R 



f du'${u',u>)jTF n {u)F n {u')^± forO<v<R, 



ip(u, V, to) 

1/4 [• ^ tt ( \ 

$(u,v,u) = ^= y du'#u» J] F n (u)F n (u') for R < v. 



(18) 



(19) 



(20) 



(21) 



III. LINEAR CASE: A = 

In the linear case (A = 0), Eq. (20) assumes the form 



oc 



Vi (R) Ui(R) i^du . 



Taking into account that the set of functions F n {u) is complete and orthonormal, one can rewrite Eq.(22) as 



^{ucj-^^UiiR) Vt(R)- 1)^ = 0, 



where 



Fi(u) <j>(u, to) du. 
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(22) 



(23) 



(24) 



Thus, we can conclude that the solution of the linear eigenvalue problem can be presented in the form 

<Pi=6 kn Af n (25) 



(26) 



Eq. (26) determines the frequency (u> n = A„) of the n-th eigenstate and Eq. (25) with Af n being a normalization 
constant and Skn being the Kronecker delta, yields its amplitude. 

Introducing Eqs. (24) and (25) into Eqs. (21) we obtain that the eigenfunction $„(u, v) = ^(u,v,u} n ), which 
corresponds to the eigenvalue given by Eq. (26), can be expressed as: 



$„(«,«) =M n F n {u) [ M^-9(R-v)+ ^M-6{v-R) 



V n {R) 



U n {R) 



(27) 
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For even values of n: n — 2m (to = 0, 1, 2, ..), Eq. (26) always has a solution and for vR — > and the eigenvalue is 
given by 

A2m ^ ( r(r l! 2) ) 4(1+^)2- (28) 

For n = 2m + 1, to = 0, 1, 2, ... the bound state exists only for vR > 1 and near the lower bound the energy of the 
bound state is given by 



to! \ {vR - l) s 
2T(m + |) / v 2 R 3 



A2m +! ~ I ^F7Z , 3^ ) ..2 pa • ( 29 ) 



In the limit of large radius of curvature R and moderate n, i.e., vR ^> n + i , we obtain from Eq. (26) that the 
eigenvalues A„ are determined by 

A„ = ^l-^V (30) 



2 V f# . 

Thus the bound state energy decreases when the curvature of the chain increases and in the limit fl^oowe obtain 
the straight-line result: A = v/2. 

It is interesting to return to the Cartesian coordinates and to consider the shape of the bound state wave function. 
Let us consider the case n = 0. It is seen from Eqs. (A3) and (A9) that 

$o = A/o e- Ao » (crfc (v%i?) - ^) + crfc (yyfc) v)) ^ ( 31 ) 

where the function v(x, y) is given by Eq. (9) When R — > oo (straight waveguide) the wave function is localized in the 
y-direction only. However, for finite R the function is localized both in x- and y-directions. The localization length 
in the x-direction is proportional to R. The expression of Eq. (31) will also be used as a starting point in our direct 
numerical simulations of Eq. (10); see section V below. 

IV. NONLINEAR CASE 

In the following, we restrict ourselves to the case of the waveguides with small or moderate curvature (l/R < 1). 
We use Darwin's expansion of the parabolic cylinder functions [30] : 



U(a,x)V(a,x)*J {x2+4ay x 2 +4a»l. (32) 



Using this approximation, Eq. (20) may be represented in the form 

Vr 



vw 1 ' i ^F{\ct>\ 2 )(j){u,u) = ^ f F n (u) F n (u')y/u 1/2 R + 2n + 1 <j>{v! » du' , 

2 n=0 J 

or taking into account that 



(33) 



in the equivalent form 



-^ + (1 + ^)^-2 cf>(u,Lu) + — \<f>\*<f>{u,u>) = 0. (34) 



Thus, eliminating the waves in the linear medium in which the waveguide is embedded and applying the inverse 
Fourier transformation with respect to Eq. (11), leads to the following equation for the waveguide function: 
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(f>(u, t) = ip(u, R, t) 



(35) 



y-dl-i(l + ^)d t -^j 0(M) + ^(M)| 2 <KM)=O. (36) 

Thus, the dynamics of the system is described by the pseudo-differential or, in other words, the nonlocal, in time and 
space, equation. The nonlocal character of the nonlinear waveguide dynamics is due to the existence of two paths 
for the excitation energy transfer: directly along the waveguide and via the linear medium in which the waveguide is 
embedded. 



A. Stationary states 

We are interested here in the stationary solutions of Eq. (36) of the form 

<j>(u,t) = $(w)e jA2 *. (37) 
Here the shape function <!>(«) satisfies the equation 



-dl + (1 + ^) X 2 \ j $ + ^|$| 2 $ = 0. (38) 

Let us consider in more detail the case of repulsive excitations (A > 0). This case is the most interesting from the 
point of view of the interplay between the nonlincarity and curvature, because, for the straight waveguide (R — > oo), 
Eq. (38) has no localized solutions. 

When the nonlinearity parameter A is sufficiently large, one can use the so-called Thomas-Fermi approximation 
[31] in which one neglects the term d 2 in Eq. (34), and one finds a density profile 

- ^ - v/T+^W) , -uo<u<u , 

\<f>(u)\ 2 =0, \u\>u , (39) 

where u = R^J — 1. 

By using Eqs. (37) and (10), it is straighforward to show that in terms of the channel wave function (18) the 
Hamiltonian (1) and the norm (4) may be written as 



oo 

\2 



i r V2 

H = --vAR / 4> 4 (u)du- — N, (40) 



N = ~uAR^ j <>H«)<h,. (41) 

Inserting Eq. (39) into Eqs (40) and (41), we obtain the following expressions for the energy and the norm of the 
nonlinear excitations: 



vR ./ 16A\ / v 2 18A . , , , , 



v a#T^\Ta7 1 ™M\Ta- ' 



From Eqs (42) and (43) one can obtain that for NvA/R > 1 



6 



H 



6 \ A J V 2 



while for NvA/R < 1 



4 20 \4RJ 

B. Non-stationary dynamics 

When R — > oo Eq. (36) assumes the form 

(V-^-zft-^) 0(u,t) + ^|0(u,t)|V(«, «) = <>, (44) 



which is the nonlinear Hilbert-NLS equation introduced in Ref. [32]. 

In the limit of weak nonlinearity (A <C 1) and small curvature (R ^> 1) one can significantly simplify Eq. (36) by 
using the Ansatz 

(f)(u, t) = exp (i t v 2 /4) ip(u, t) (45) 
and assuming that <p(u,t) depends slowly on u and t. Inserting Eq. (45) into Eq. (36), we obtain 

^-^+(l + ^)(^-ia t )-^ p(«,t) + ^( u ,i)|V«, *) = <>. (46) 
Considering the scaling 

t = te~ 2 , u = ue~ 1 , R = R€~ 2 , A = Ae 2 (47) 
for e — > 0, Eq. (46) reduces to 

8 2 v 2 u 2 v 2 A^ l2 

l df = -M 2 - iP+ Jt iP+ — Mip - (48) 

Thus the behaviour of nonlinear excitations in a curved waveguide is equivalent to the behaviour of nonlinear exci- 
tations in the parabolic potential whose curvature coincides with the curvature of the waveguide. In the linear case 
A = the corresponding eigenvalue problem gives the eigenvalues presented in Eq. (30). Furthermore, it is interesting 
that this reduction gives rise to an effective one-dimensional Gross-Pitaevskii (GP) equation, analogous to the one 
used in the study of Bose-Einstein condensates in cigar-shaped traps [33] . While we do not pursue this analogy in 
detail in the present proof-of-principle paper, we note that it would naturally be of interest to examine excitations 
known in the context of the GP equation, such as e.g., dark solitons [34] and their dynamics, in the present setting. 

V. NUMERICAL RESULTS 

We start by demonstrating the results of the linear case of Eqs. (26) and (31). The eigenvalue (energy) of the linear 
case as a function of curvature is shown in Fig. 1, while the lowest energy, bound state wavefunction of the linear 
problem is given in Fig. 2. 
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FIG. 2. The wave function of the lowest bound state of the linear problem (^4 = 0) is shown in the figure for R = 10, v = 1 
(obtained through Eq. (31)). 

In order to demonstrate that this linear bound state persists in the nonlinear limit we have performed full dynamical 
evolution simulations of Eq. (7), with an initial condition of the form expressed in Eq. (31), and demonstrated in Fig. 
2. We note in passing that similar results have been obtained with other initial conditions such as e.g., a Thomas- 
Fermi initial profile of the form of Eq. (39). In particular, we show typical numerical simulation results in Fig. 3 for 
R = 10, v = A = 1. Notice that the 8 function was represented as 




with e = 0.05. The contour plot of Fig. 3 shows the result after a numerical evolution of 100 time units of Eq. (7). 
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FIG. 3. Contour plot of the solution of the partial differential equation of Eq. (7) with f(x) — x 2 /(2R), R = 10, v = A = 1 
and initial condition given by the linear profile of Eq. (31). 



The dynamical development indicates that after an initial transient the original linear profile slightly reshapes itself 
into the nonlinear solution depicted in Fig. 3. In the process, some radiation waves ("phonons") are shed, that are 
absorbed by the absorbing boundary conditions used in a layer close to the the end of the domain (our computational 
box is of size 25 x 25). 

Beyond the proof-of-principle simulations for various initial conditions theoretically derived in sections III and IV, 
we also attempted to examine the dynamics of the nonlinear excitations of the channel. This was done using the 
following numerical protocol: after obtaining a quasi-relaxed nonlinear localized mode for the channel of the form 
y = x 2 /(2R), we moved the channel to a new position, namely y = (x — l) 2 /(2R). Notice that similar in spirit 
experiments have recently been carried in Bose-Einstein condensates [35], where the magnetic trap confining the 
condensate is displaced to a new position and the ensuing dynamics of the condensate are observed. The position 
of the center of mass of the initial condition profile was approximately obtained (using a trapezoidal approximation 
to the relevant two-dimensional integrals) as (x(t — Q),y(t — 0)) = (0.178,0.862). The new bottom of the channel 
(hence the point to which the center of mass should approach) in this case is (x,y) = (1,0). In Fig. 4, the process 
of relaxation to this new equilibrium is shown as a function of time for a very long dynamical simulation of t up to 
1000 time units. In this run, we observe (after an initial transient) a slow relaxation towards the new minimum of 
the potential well. Notice that despite the Hamiltonian nature of the model, the excitation of an "internal mode" of 
the nonlinear wave [36] can be dissipated due to mechanisms of coupling to extended wave, phonon modes, such as 
the ones reported in [37]. 
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FIG. 4. The top panel shows, for times between t = 500 and t = 1000, the time evolution of the center of mass in the x — y 
plane (the motion proceeds from top left to bottom right as time evolves) . The middle panel shows the time evolution of the 
x-position of the center of mass of the nonlinear excitation, while the bottom panel the slow evolution of the y-position of the 
center of mass. 

VI. SUMMARY 

In summary, we have shown that: 

• In two-dimensional media with a curved infinitely thin waveguide (quantum channel) there exist bound states 
for linear and nonlinear sclf-intcracting excitations; 
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• The finite curvature of the waveguide provides a stabilizing effect on otherwise unstable localized states of 
repelling excitations; 

• The binding energy of both linear and nonlinear localized excitations decreases when the curvature of the 
waveguide increases. 

• Such linear bound states as the ones found here persist in the nonlinear dynamical problem as localized exci- 
tations. These have been found to be robust for different initial conditions and are centered at the minimum 
(point of largest curvature) of the parabola. When, the mode is initialized away from this minimum, it slowly 
relaxes to it. 

While our motivation was originally provided by the embedding of a waveguide in a two-dimensional photonic crystal 
(as well as from more general geometric considerations), interesting analogies have arisen through our investigation, 
that warrant further study. In particular, we note the analogy in the weakly nonlinear regime of the equation for the 
waveguide /quantum channel with that of the Bose-Einstein condensate behavior in the presence of a magnetic trap 
in atomic physics. Another topic worthy of further investigation is a more detailed numerical study of the stability 
of the nonlinear localized modes we have identified. Finally, another direction that could be of further interest is the 
examination of a case of a finite (rather than infinitesimal) width channel [which can be computational achieved e.g., 
by allowing the parameter e of Eq. (49) to vary towards larger values]. The examination of thresholds for genuinely 
two-dimensional instabilities, such as e.g., the transverse or the snaking instability, would be of particular interest 
within the latter context. 

Such studies are currently in progress and will be reported in future publications. 
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APPENDIX A 



To find an expression for the Green's function G(u, u'; v, v') we expand this function in terms of of the eigenfunctions 
of the equation 

0-^ = -„F, (Al) 

under the boundary conditions 

F^O, u^±oo. (A2) 

Eq. (Al) under the boundary conditions (A2) has the solutions 

F n (u) = a n e-° 2 / 2 H n (a), n = 0,1,2,... (A3) 

((7 = u yJ\/R), which correspond to the eigenvalues 

^ = 2n+l. (A4) 

Here H n (z) is the Hermite polynomial [30] and a n = (X/Rir) 1 ^ 4 (2™ n!)~ 1,/2 is the normalization constant. In terms 
of the functions (A3), the Green's function has the form 

oo 

G(u, u'; v, v') = F n (u) F n {u') g n (v, v'). (A5) 

n=0 

Inserting Eq. (A5) into Eq. (14) for the coefficients g n (v,v') we obtain the equation 
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d 2 g n \ 2 A 



(lr 2 R2 v 2 gn~^(2n+l)g n = 6(v-v') 7 (A6) 
under the condition 

g n -> 0, t; -> oo. (A7) 

The line v = 0, which is the y > R/2 of the ordinate axis, is used as the cut [38]. Following the usual machinery (see 
e.g. [38]), to avoid discontinuities in the solution along the cut we shall require that 

= at v — when n = 2m, 

civ 

g n (0) = when n = 2m + 1 (m = 0, 1, 2, ..). (A8) 
The solution of the equation (A6) with regard to the boundary conditions (A7) and (A8) has the form 



^ {V(n + 1 v VWR) U(n + l,v> ^/R) 9(v> - v) + 



U(n+^,v^2X/R)V(n+^v^2X/R)0( v - v')} (A9) 



where U(a, x) and V(a, x) are the Weber parabolic cylinder functions [30] and 0(v) is the Heaviside step-function. 
Another equivalent representaion in terms of Whittaker functions [30] is 



1 . 2 -™/ 2+3 / 4 n! , / n 111.,. 
y(n+ 2'^ = ^r (i + l) M (-2-4'4'r 
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